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Abstract 

QCD at a finite quark-number chemical potential /z has a complex fermion determinant, which 
precludes its study by standard lattice QCD simulations. We therefore simulate lattice QCD at fi- 
nite fj, in the phase-quenched approximation, replacing the fermion determinant with its magnitude. 
(The phase-quenched approximation can be considered as simulating at finite isospin chemical po- 
tential 2/z for Nf/2 n-type and Nf/2 d-type quark flavours.) These simulations are used to study 
the finite temperature transition for small fi, where there is some evidence that the position (and 
possibly the nature) of this transition are unchanged by this approximation. We look for the ex- 
pected critical endpoint for 3-flavour QCD. Here, it had been argued that the critical point at zero 
fi would become the critical endpoint at small ^, for quark masses just above the critical mass. 
Our simulations indicate that this does not happen, and there is no such critical endpoint for small 
fi. We discuss how we might adapt techniques used for imaginary fi to improve the signal/noise 
ratio and strengthen our conclusions, using results from relatively low statistics studies. 
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I. INTRODUCTION 



Relativistic heavy-ion colliders allow one to study hadronic and nuclear matter at high 
temperatures where it undergoes a transition to a quark-gluon plasma. While the highest 
energy colliders (RHIC and the forthcoming heavy-ion program at the LHC) study only the 
very low density regime where the baryon-number density is too small to have much effect 
on the thermodynamics, lower energy relativistic heavy-ion colliders can probe the region 
where baryon-number density is appreciable. 

For physical u, d and s quark masses, the finite temperature transition at zero baryon- 
number density is predicted to be a rapid crossover rather than a true phase transition 

m 

. It IS expected that, at high enough baryon-number densities, this transition will 
become first order. The point at which the change from a crossover to a first-order transition 
occurs would be a critical point, expected to be in the universality class of the 3-dimensional 
Ising model. This critical point is referred to as a critical endpoint, and is expected to be 
the most interesting feature of this intermediate density regime of the QCD phase diagram. 

While finite temperature QCD is straightforward (but tedious) to simulate on the lattice, 
QCD at a finite quark-number chemical potential /i has proved intractable. The reason 
is that at finite /i the fermion determinant becomes complex, with a real part having an 
indefinite sign. Since all the standard lattice QCD simulation methods rely on importance 
sampling, they fail for such systems. 

In the region of small /i, close to the finite-temperature phase transition, methods have 
been developed to circumvent this sign problem. These methods fall into several classes. One 
such method involves simulating lattice QCD at a carefully selected set of parameters where 
no such sign problem exists and using the ratios of determinants to reweight to the region 
of interest 4|- Such multiparameter reweighting only works provided there is significant 
overlap between those configurations which are important for the chosen set of parameters, 
and those which are important for the original set of parameters. A second class of methods 
are those which rely on analyticity in /i or related parameters. These include series expansion 
methods 5|, |6| , which expand the Boltzmann weight and the observables as power series in 
/i, calculating the coefficients in simulations at zero /i. Since the higher order coefficients 
require the calculation of higher order fluctuation quantities, this ultimately limits their 
utility. Other analyticity methods involve simulating in a domain of parameters such as at 
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imaginary /i, where there is no sign problem, and analytically continuing the results to the 
desired domain (in this case, real yu), typically by fitting the results to a power series [3, ^■ 
There exist variants where different parameters are used for the analytic continuation such as 
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9j. Another way of avoiding the sign problem is to use canonical methods [10|, lUl, ll2[. Here 
the sign problem is encountered in Fourier transforming to obtain the canonical ensembles 
at fixed quark number. 

We have adopted the alternative approach of ignoring the phase of the determinant and 
replacing the determinant by its magnitude. This can be thought of as simulating QCD with 
Nf/2 u type quarks and Nf/2 d type quarks, with a chemical potential /i/ = 2/i for isospin 
(J3). {Nf is the number of quark flavours.) For low temperatures, there is a critical point 
III = fic above which the system enters a superfiu id p hase, with a charged pion condensate 



which breaks I3 symmetry spontaneously [l3 



14 



15|. At zero temperature fic = 'fn-w- Since 



this phase does not exist for QCD at finite /x, the phase-quenched approximation breaks 
down at the boundary of this superfiuid domain, if not before. 

nn 

The Taylor series calculations of the Bielefeld-Swansea collaboration p, 116( showed ev- 
idence that the /i dependence of the transition temperature Tc for full 2-fiavour QCD was 
similar if not identical to that of the phase-quenched approximation (finite isospin chemical 
potential), at small /i. The /x dependence of this transition for full 2-fiavour QCD, ob- 
tained from the imaginary quark-number chemical potential simulations of de Forcrand and 
Philipsen, [3] was consistent with being identical that observed in our direct simulation of 
the phase-quenched theory A random matrix model of 2-fiavour QCD at finite temper- 
ature and chemical potentials also predicts that the dependence of on quark-number and 
isospin chemical potentials should be identical for /i < 771,^/2 In addition Nambu-Jona- 
Lasinio models for QCD have transition temperatures which exhibit the same dependence 
on quark-number and isospin chemical potentials for /i < 7X1^^/2 [l^, 2^. This strongly sug- 
gests that the /i dependence of Tc is the same for phase-quenched and full QCD for small 
/i. We shall indicate later that the simulations discussed in this paper are consistent with 
this assumption. There is, however, one lattice result which contradicts this assumption. 
The 3-fiavour calculations of the Bielefeld- Swansea collaboration indicate that while the 
dependence on quark-number and isospin chemical potentials are consistent at larger quark 



masses, they are not at small quark masses 
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22l |. However, because of large statistical 



errors, the observed difference in slopes was less than two standard deviations. In addition. 
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these simulations were performed using the R algorithm, which could potentially introduce 
larger than expected updating errors, due to the discretization of molecular-dynamics 'time', 
in the fluctuation quantities used to obtain these results. 

For 3-flavour QCD at zero chemical potentials, the finite temperature transition is first 
order at small quark mass m. For larger m the transition softens to a crossover with no phase 
transition. At m = rric, where the nature of the transition changes, the finite-temperature 
transition is a critical point in the universality class of the 3-dimensional Ising model 
Similar behaviour is seen for 2 + 1-flavour QCD, and for the physical u, d and s quark 
masses, the transition is predicted to be a crossover j^, 3|. It has been suggested that rric 
would increase with increasing /i, becoming the critical endpoint. If so it should be possible 
to tune this endpoint to be as close to yU = as desired by choosing m just above rric- 

Hence we simulate 3-fiavour lattice QCD at /i/ < m^, for several masses close to ttIc, and 
determine the nature of the finite-temperature phase transition using fourth-order Binder 
cumulants for the chiral condensate. For these studies we use simulations on 8'^ x 4, 12^ x 4 
and 16^x4 lattices. Our simulations indicate that there is no critical endpoint for m > mc(0), 
and mc{^i) actually decreases (slowly) with increasing /i. Preliminary results from these 
simulations have been reported at various conferences, the most recent being Lattice2007 
23|. This absence of the expected critical endpoint at small has been observed by de 



Forcrand and Philipsen using analytic continuation from imaginary [2|. Our simulations 
use the exact RHMC algorithm 2J], since in the inexact hybrid molecular-dynamics used 
in earlier simulations, the Binder cumulant had such strong dt^ dependence as to lead to 
incorrect conclusions as to the nature of the transition 251 ]. 

The relatively weak dependence of the Binder cumulant on /ij and the statistical errors 
in determining it, even in high statistics runs, make it difficult to determine the sign of 
the slope dB^^/djji'] and hence dmc/dfij with certainty. Similar difficulties occurred with the 
methods of de Forcrand and Philipsen, but they were able to calculate the slope directly with 
much higher precision, using reweighting methods 26|]. We have performed studies which 
indicate that similar methods should work for the phase-quenched simulations. However, 
on the larger lattices we use, it is unclear whether these methods will be significantly more 
efficient than simply increasing statistics. As yet, we have insufficient statistics to achieve 
results for the slope of the Binder cumulant. However, we are already able to determine the 
slope of Pc- 
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Section 2 describes phase-quenched lattice QCD. In section 3 we describe our simula- 
tions and results. Exploratory studies of reweighting techniques are described in section 4. 
Section 5 is devoted to discussions and conclusions. 



II. PHASE-QUENCHED LATTICE QCD 

Phase-quenched lattice QCD with eight staggered quark flavours (or two staggered quark 
fields, each with four 'tastes') has the fermion action 

sites 

where Jp^^r^fij) is the standard staggered quark transcription of ^with the links in the +t 
direction multiplied by exp(|r3yU/) and those in the —t direction multiplied by exp(— ^rs/i/). 
Since we are performing simulations outside of the superfluid phase, the explicit symmetry- 
breaking interaction of our earlier studies is unnecessary. 

To simulate Nf flavours using the RHMC algorithm, this is replaced by the pseudo- 
fermion action 

S,f=plM-''f/% (2) 



where are the momenta conjugate to the pseudo- fermion field '?/'[36|. Here 



M = [pi^fii)+m]^[pi^f,i) + m]. (3) 
In the RHMC algorithm A4~^f/^ (and Ai'^^f^^^) are replaced by rational approximations. 



using a speculative lower bound 25|. It is interesting to note that these rational approxi- 
mations provide similar infrared protection to what a symmetry-breaking interaction would 
give. 

For 8 flavours, and indeed for any even number of flavours, fij has the interpretation of an 
isospin chemical potential, for a theory with Nf/2 -u-type quarks and Nf/2 d-tjpe quarks. 
Since we are interested in this phase-quenched theory as an approximation to QCD with a 
quark- number chemical potential fi = fii/2, we are free to choose any integral Nf. In fact 
we shall work with Nf = 3. 

As we have shown in earlier work, the Binder cumulant which is used to extract the 
nature of the finite-temperature transition is very sensitive to the updating increment dt 
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in the older, inexact, hybrid molecular-dynamics (R) algorithm 25|. This is the principal 
reason that we have switched to the RHMC algorithm. 

As mentioned in the introduction, such theories are known to undergo a phase transition 
to a superfluid phase with a charged pion condensate and orthogonal charged pion excitations 
which are true Goldstone bosons at low temperatures, as fij is increased. At zero temperature 
this transition occurs at fii = fic = fn-K- At high enough temperatures the system should be 
in the quark-gluon phase for all /i/, and no such transition is expected. 



III. SIMULATIONS AND RESULTS 

We perform simulations of 3-fiavour lattice QCD at finite fij and temperature on 8^ x 4, 
12^ X 4, and 16^ x 4 lattices. We use rational approximations to Ai^^^^ and A^^^/^^ in 
these RHMC simulations which are valid provided the spectrum of A4 is in the range [1 x 
10~^, 25]. (2 runs were made using smaller speculative lower bounds for testing purposes.) 
We performed runs with quark masses m = 0.02, m = 0.025, m = 0.03 and m = 0.035 on 
8'^ X 4 and 12^ x 4 lattices. At the lowest mass, we only ran simulations for fij = 0. For the 
other 3 masses we ran simulations at fij = 0, fij = 0.2 and fij = 0.3. In addition, we ran 
simulations on 16^ x 4 lattices at m = 0.03 at all 3 /i/ s and at m = 0.025 with fij = 0. The 
masses are chosen such that the lower 2 masses lie below rric and the higher 2 masses lie 
above rric. The choice of /ij values is to cover the region < fij < m^, where is estimated 
to lie in the range 0.4 ^ m^r < 0.5 for 0.025 < m < 0.035. (This comes from measurements 
of for m = 0.03, fii = at j3 = 5.10 and (3 = 5.15, which bracket the region of interest, 
on 8'^ X 16 and 12^ x 24 lattices. Estimates for other m values were made using PCAC.) 

For our 12^ x 4 simulations, where we have the highest statistics, we have run for 300,000 
length-1 trajectories for each of 4 (or more) (3 values close enough to the transition to access 
this transition using Ferrenberg-Swendsen reweighting in j3, at each {rrijfii). For the 8^ x 4 
and 16^ x 4 simulations we have performed 300,000 trajectory runs at each of 2 /5s at each 
(m,/i/). We have made 5 independent stochastic estimates of the chiral condensate ipil) 
and the isospin density jg = dSf/dfij after each trajectory, to enable us to make unbiased 
estimates of the susceptibilities and Binder cumulants. 

For any observable X, the susceptibility xx is defined by 

Xx = ^(x'-{Xf), (4) 
6 



where V is the spatial volume and T = 1/Nt is the temperature. The overlining of X 
indicates that these are lattice averaged quantities. The fourth-order Binder cumulant for 
X is defined by 

- {(x-(x)W 



27l |. These quantities are measured at the value P = Po of the simulation and extrapolated 



to nearby /3s, by Ferrenberg-Swendsen reweighting: 



(exp[-6(r/T)(/?-/jo)^o]X)^, 
^ {exp[-6{V/T){(3-f3o)Su])p, ^ ^ 



28l | where 

Sa = l- ^Re TiaUUUU (7) 

appropriately averaged over the lattice and over plaquette orientations. The position of the 
transition, (3c, can be estimated as that of the peak of the susceptibilities, or the minimum 
of the Binder cumulants. We have noticed that the measured jScS from the susceptibilities 
of various observables and from the Binder cumulants are close. 

The Binder cumulant for the chiral condensate is used to probe the nature (as well as 
the position) of the transition. In the infinite volume limit, S4 = 3 at a crossover, 3^ = 1 
at a first-order transition and B4 = 1.604(1) at a 3-dimensional Ising critical point. If there 
indeed were a critical endpoint, for m > mc{0) B4 would start at a value above the Ising 
value for /i/ = (close to 3 for really large lattices) and decrease, passing through a value 
close to the Ising value at the critical endpoint, eventually approaching 1 for large fij. On 
large enough lattices, finite-size scaling predicts that lines of versus /i/ for different lattice 
sizes will cross at the Ising value. Similarly lines of B4 versus m for different size lattices 
will cross at the Ising value as m is varied. 

In figure [1] we plot the Binder cumulants at the transition point as functions of /if for 
m = 0.025, m = 0.03 and m = 0.035 for the various lattice sizes. For m = 0.035 B4 on 
the 12^ X 4 lattice starts at a value significantly above that for an Ising critical point and 
appears to increase with increasing fij, and hence shows no evidence for a critical endpoint 
the slope of this straight line fit is 0.68(26). Similarly, for m = 0.03, BA starts above the 
Ising value and increases with increasing fij on the 12'^ x 4 and 16^ x 4 lattices. For the 
12^ X 4 lattice the slope is 0.39(22), and for the 16^ x 4 lattice, this slope is 0.76(53). For 
m = 0.025, there is no evidence for any fi'j dependence for B4 on the 12^ x 4 lattice and 
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it remains below the Ising value for the range of fij considered. We note that the 8^ x 4 
slopes appear negative for m = 0.035 and m = 0.025 and positive for m = 0.03, which we 
interpret as meaning that we have insufficient statistics to determine the very small slopes 
of the 8^ X 4 lines. 
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FIG. 1: Graphs showing the ^'j dependence of the Binder cumulants for the chiral condensate tpip 
at the transition: a) for m = 0.035, b) for m = 0.03, c) for m = 0.025. The dashed line is at the 
Ising value. 

None of the slopes we have measured is much more than 2| standard deviations from 
zero. However, the fact that the two 12^ x 4 slopes and the one 16^ x 4 slope for m > rric 
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are all positive makes it less likely that this is a statistical fluctuation. For m < rric we can 
draw no conclusions. 

We now turn our attention to the mass dependence of -B4 at fixed fij values. Figure [2] 



shows the m dependence of i?4 for /i/ = 0, /i/ 
SU(3) Nf = 3 ycxi = A=0 



0.2 and fij = 0.3. First we note that 
SU(3) Nf = 3 /Xi = 0.2 A = 
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FIG. 2: Graptis showing tfie m dependence of the Binder cumulants for the chiral condensate il^ijj 
at the transition: a) for /i/ = 0, b) for /x/ = 0.2, c) for fij = 0.3. The dashed line is at the Ising 
value. 

the intersection of the curves for the different lattice sizes intersect at B4 close to its value 
for the 3-dimensional Ising model. This is strong evidence that this critical point is in 
the universality class of the 3-dimensional Ising model, as predicted. We therefore use the 
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masses for which the 12^ x 4 Binder cumulants achieve the Ising value as our estimate for 
the position of the critical point for the /i/ under consideration. We get mc(0) = 0.0265(3), 
mc(0.2) = 0.0259(5) and 771^0.3) = 0.0256(4). A straight line fit yields 



This suggests that rric decreases with increasing rather than increasing as would be needed 
for a critical endpoint. Note also that if we were to use the intersections of the curves for 
different lattice sizes as our estimates for rric, this would slightly increase our estimate for 
mc(0), slightly decrease our estimate of mc(0.3) and leave our estimate of mc(0.2) essentially 
unchanged. This would make the slope even more negative. In addition, since the transition 
temperature decreases with increasing rric in physical units will decrease slightly faster 
than the rric in lattice units, which we have presented here. 

We have also examined the Binder cumulants for the isospin density j'q, and find that 
these are consistent with those for the corresponding chiral condensates. However, since our 
estimates for j'q are much noisier, the errors on -B4(jo) are considerably larger than those 
for i?4 ('?/;'?/'), which makes them less useful. The Binder cumulants for the plaquettes are 
appreciably larger, which is expected, since these should be a reasonable approximation 
to the energy-like order parameter whose Binder cumulant would approach 3, even in the 
first-order regime and at the critical point. 

The critical behaviour of this theory will be described by an effective Hamiltonian which 
is a linear combination of 3 fields, each of which has finite size scaling properties with 
one critical exponent. One such field has the scaling behaviour of a 'magnetization', a 
second that of an 'energy' and the third that of a 'density'. Each will be given as a linear 
combination of ipip, Sg = {1 — ^TinUUUU) and Jq. The sirnpler case at /i/ = 0, where there 
are only 2 fields to consider, has been studied in reference pj. In that paper they were able 
to find simple expressions for the two 'mixing' parameters. In our case there are six such 
'mixing' parameters, we have been unable to find the six equations required to determine 
these coefficients. If we were able to obtain these eigenmodes of the renormalization group, 
the Binder cumulant of the 'magnetization' order parameter would pass through the Ising 
value, once our lattice is large enough that subdominant terms in finite size scaling relation 
could be ignored. Until we can find such relationships we use the fact that, on large enough 
lattices, the magnetization component of the chiral condensate will dominate and its Binder 




(8) 



10 



cumulant will approach that of this eigenmode. The fact that the Binder cumulants for 
the chiral condensate cross close to the Ising value is evidence that this field is not strongly 
affected by contamination from non-magnetic eigenmodes, on the lattice sizes we use. 

Using Ferrenberg-Swendsen reweighting again, we calculate the chiral susceptibilities and 
measure the positions and values of the peaks. We observe that the positions of these peaks 
are very close to the minima of the Binder cumulants. Finite size scaling tells us that, at 
the critical point. 



where L is the spatial extent of the lattice and is the critical temperature. Hence if we 
plot L~ixiP'ii){L,Tc) as functions of m for different L values, the curves should cross at the 
critical point. In figure [3] we plot this quantity for L = 8 and L = 12, for each of our fij 
values. Here we have taken 7 = 1.237 and u = 0.630 as the required critical indices for the 
3-dimensional Ising model. 

Because it is clear that the points on this graph do not fall on straight lines and the 
curves for different lattice sizes cross at rather shallow angles, a quantitative estimate for 
the position of the crossing would be difficult to obtain. What is clear is that the curves 
for the different lattice sizes cross somewhere between m = 0.25 and m = 0.3 for fij = 
and fii = 0.2 and close to m = 0.25 for fij = 0.3, which is consistent with our estimates of 
rricifii) from Binder cumulants. 

As well as trying to determine the nature of the finite temperature transition as a function 
of /i/, and measuring observables and susceptibilities, the positions of the minima in the 
Binder cumulants, and the positions of the maxima in the various susceptibilities yield 
predictions for /3c the transition (3 values. The /x/ dependence of Pc will ultimately yield the 
/i/ dependence of the transition temperature Tc. This not only requires that we know Tc at 
Hi = 0, which we can obtain from the numerous measurements by other groups, but it also 
requires that we know the renormalization group running of (3 with lattice spacing a. On 
the coarse lattices we use, 2-loop perturbative running of the coupling constant which has 
been used earher, is clearly suspect. Hence we present only the /x/ dependence of Pc in this 
paper. Associated with our present simulations aimed at determining the equation-of-state 
for phased-quenched (lattice) QCD, we will measure the running of jS directly with the same 
action and masses as are used here, on zero temperature lattices. At that time we will be 
able to predict the fij dependence of T^.. 
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yU,i = chiral susceptibility 



yu,i = 0.2 chiral susceptibility 
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FIG. 3: The rescaled chiral susceptibilities x,/;^ a-s functions of m: a) for jij = 0, b) for fij = 0.2, 
c) for fij = 0.3. 

In figure H] we plot the measured values of jSc against fi'j for each of the quark masses. 
Straight line fits appear adequate with our current statistics. Although better fits could be 
obtained with a fij term for m = 0.025 and m = 0.035 - the m = 0.03 straight line fit is 
excellent - the coefficients are clearly very small, and with only 3 points on each curve, such 
an exact fit is hard to justify. These fits are to the more extensive 12^ x 4 'data'. We have 
plotted the 16^ x 4 points on the same graph. These indicate that the finite size effects on 

are very small. For comparison with the work of others, these fits are: 

/3e = 5.13418(10) - 0.1743(18)/i^ m = 0.025 (10) 
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/3c = 5.14385(8) - 0.1711(13)/i? m = 0.030 (11) 
/3e = 5.15326(10) - 0.1735(16)/x? m = 0.035 (12) 

and f3c = 5.12377(10) at m = 0.02, fn = 0. 

SU(3) Nf = 3 A=0, jSs at minima of B4 
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FIG. 4: Transition /3, /3c as functions of /i| for chosen masses. The hnes from top to bottom are 
for m = 0.035, m = 0.03 and m = 0.025. The isolated point is for m = 0.02. 

IV. REWEIGHTING STUDIES 

As we saw in the previous section, the weak dependence of the Binder cumulants on 
/i/, and the sizable statistical errors in determining this fluctuation quantity mean that the 
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observation that B4 increases with /ij, while strongly suggested is not definitive. Similar 
difficulties arise for simulations at imaginary fi. Here, de Forcrand, Kim and Philipsen have 



circumvented this difficulty by calculating dB^/dfi^ directly 26|. They do this by calculating 
B4^{fi) and B4^{fi + 6fi) in the same simulation. This is achieved by including the ratio of 
determinants 

p = det[M{fi + 6pff^'']/det[M{ijff^^] (13) 

as a weight in the measurement of tpiplfi + Sn) from the same ensemble at chemical potential 
II as is used to measure 'ijj'ijj{fi), namely 

w-v<A.+»)w.= <'"^*'^+'''»^ (14) 

and similar expressions for higher powers of ipip{fi + Sfi). 

Since exact calculation of such determinants is expensive, de Forcrand et al. used unbiased 
stochastic estimators for the ratio of determinants, in particular, 

p = (exp[-r7tA^(;,)-^//i6^(^ + + 5p)''f'^^M{p)-''''^S + (15) 

where r] is Gaussian noise. The advantage of this method is that any finite number of noise 
vectors gives an unbiased estimator of the determinant. These authors reweighted from 
= and performed a multistep reweighting to = 0.1 

Whereas it appears that de Forcrand et al. limited themselves to 8^ x 4 lattices, we 
are investigating applying this to 12^ x 4 lattices, since the results of the previous section 
make it unclear whether the slope OB^^/dp^i is the same for 8'^ x 4 lattices as it is for larger 
lattices. We first investigated the possibility of reweighting from /// = to /x/ = 0.1 in a 
single reweighting, but analysis of a few configurations quickly convinced us that although 
the overlap might be reasonable, the fiuctuations were so large as to make it impossible to 
obtain a reasonable estimate of the determinant without use of far more noise vectors than 
is reasonable. We then went back to a reweighting from pj = to pj = 0.01 as a single 
step process, and one that could be used as a basis for a multistep reweighting to an even 
larger pj. For this trial run we used 1500 configurations at m = 0.03, separated by 200 
trajectories. For each configuration we used 200 noise vectors with 6pi = 0.01 and the same 
set of noise vectors with 6pi = —0.01, making use of the fact that the determinant for a 
single configuration remains unchanged under pi —* —pj to remove 0{6pj) fiuctuations in 
our noisy estimator. We used 1000 noise vectors for our noisy estimators for and 
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the same set for (0.01) and 'tpipl—O.Ol). This effectively removes the errors in using 
noisy estimators of the condensate from consideration. The resulting estimate of dB^/dfij 
is 3.1 ± 4.1 compared with the estimate 0.38 ± 0.22 obtained in the previous section. This 
indicates that, as expected, we need to use a S{fi'j) much greater than the 0.0001 used here, 
which will require a multistep reweighting in order to avoid large fluctuations. Our estimate 
for dPc/dfij is —0.177(9), in agreement with —0.171(1) obtained in the previous section. 
Figure [5] shows the noisy estimators of the determinants with errors that we obtained. We 
see that the errors are comparable with the difference of these determinant ratios from one 
and from their mean, which is one reason why the signal/noise ratio is so poor. 

We notice with the reweighting from fij = 0, that one problem is that the signal is of order 
5(/if) = {6fii)^, while the noise is of order 6fii. While this can be overcome with a multistep 
(multiple nis) reweighting, an alternative way of avoiding this difficulty is to start at non- 
zero /i/ where for small 5/1/, 5/1/ and 5(/i/) are of the same order of magnitude. We have thus 
tried 1-step reweighting from 1500 configurations at fij = 0.2, with = 0.01 and hence 
with S{fij) = 0.0041. For this test we ran first with 200 noise vectors for each configuration, 
and later with 1000 noise vectors for each configuration. 1000 noise vectors were used in 
estimating the chiral condensate. Using 200 noise vectors to estimate the determinant ratio 
we obtained dB^/dnj = —0.39 ± 0.56 and for 1000 noise vectors —0.54 ± 0.45. Although 
this indicates that we still do not have enough statistics, we would only need to reduce 
the statistical errors by an order of magnitude to make a definitive prediction. d(3c/dfij 
measured in the same calculations is —0.175(2) compared with —0.171(1) calculated in the 
previous section. In figure El we show our estimates of the ratio of fermion determinants. 
Even with 200 noise vectors/configuration, the ratio of determinants is well determined. 
The statistical errors are considerably smaller than the ratio's departure from unity and, 
more importantly, considerably smaller than the range of values taken by this ratio over the 
ensemble of configurations. This presumably is why little improvement in the estimate of 
dB^/dfij is obtained by increasing the number of noise vectors from 200 to 1000. Comparison 
of the estimates of the determinant ratios for 200 and 1000 noise vectors makes us confident 
that our noisy estimates are reliable. 

In case the main problem was overlap, we reduced our 6fii to 0.005. The results, however, 
were similar. djSc/dfij was well determined, while the errors in dB^/dfij exceeded the signal. 
The determinant ratio was well determined by 200 noise vectors, the errors being much 
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FIG. 5: a) Stochastic estimates of the ratio of fermion determinants at fij = 0.01 and /U/ = on a 
12^ X 4 lattice at m = 0.03, fij = 0, P = 5.143. b) Section of graph (a) showing detail. 
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FIG. 6: a) Stochastic estimates of the ratio of fermion determinants at fij = 0.21 and fij = 0.2 
on a 12'^ X 4 lattice at m = 0.03, m = 0.2, /3 = 5.137. The crosses(red onhne) are for 200 noise 
vectors; The circles(blue onhne) are for 1000 noise vectors, b) Section of graph (a) showing detail. 
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smaller than the fluctuations in the value of this ratio from configuration to configuration. 

To get an accurate estimate of dB^/dfij will require the analysis of many more configu- 
rations. These configurations should be separated by enough trajectories as to make them 
reasonably independent or they are unlikely to improve the errors, just as increasing the 
number of noise vectors for estimating each ratio from 200 to 1000 did not significantly 
improve our errors. In fact, using 1 noise- vector to estimate the determinant ratio at the 
end of each of 300,000 consecutive trajectories at /i/ = 0.2, gave similar accuracy to using 
200 or 1000 noise-vectors for each of our 1500 configurations spaced by 200 trajectories. In 
addition, producing a single trajectory takes much less computer time than a measurement 
with (say) a 100 noise-vector estimate of the determinant ratio and a 100 noise-vector esti- 
mate of tptp. In addition we need to check whether a single- or a multi-step estimation of 
the determinant ratio is more efficient, even at fij = 0.2 where it is not forced on us by other 
considerations. 



V. DISCUSSIONS AND CONCLUSIONS 



We have studied the finite temperature transition for 3-fiavour lattice QCD with a finite 
chemical potential /i in the phase-quenched approximation, where the phase of the fermion 
determinant is set to zero, using RHMC simulations. This can be considered as studying 
lattice QCD with 3/2 up-type quarks and 3/2 down-type quarks at a chemical potential 
/i/ = 2/i for isospin (I^). As we have indicated in the introduction, in the small ^{^i) regime 
- /i/ < - there are indications that the dependence of the transition temperature on /j, 
for the phase- quenched model is similar if not identical to that for full QCD. Within the 
limitations of our statistics, which only allow us to include terms linear in fi^, we find that the 
coefficient of /i^ in the fit to j3c{^i'^)-, is within 10% and probably within 5% of that obtained 
by de Forcrand and Philipsen for the full theory by continuation from imaginary {2!, . 
This gives further evidence that the /i^ dependence of T^. is the same in phase-quenched and 
full QCD. If this is true, it is reasonable to expect that the nature of the transition will be 
the same in both theories. 

Let us briefiy digress to discuss other recent work which could have relevance to the 
connection between phase-quenched and full QCD. Some recent work of Fodor, Katz and 
Schmidt, which employs the density-of-states method, uses the phase-quenched theory as a 
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starting point for their factorized reweighting 29j. This shows a small but finite shift in (3c in 
reweighting from phase- quenched to full QCD. However, all the values which they consider 
are larger than m-^jl^ and so in the region where the two theories are no longer expected to 
be similar. There has also been extensive work on QCD at finite chemical potentials using 
a random matrix/chiral perturbation theory approach 30, l3l|]- This has indicated that the 
phase of the fermion determinant becomes much worse behaved at for /i > Although 
this work does not (yet) explain why the full and phase-quenched QCD behave similarly, it 
does indicate that the same pion modes describe the physics of each, and suggests model 
calculations which might clarify the situation. 

It was expected that the critical point at zero chemical potential would move to higher 
mass at finite chemical potential. If so, for quark masses just above the critical mass at 
fi = fii = 0, this would become the sought-after critical endpoint where the crossover 
at /i = would change to a first-order transition. Our simulations for m close to mc{0) 
indicate that this does not happen, but rather mc{^j) decreases with increasing /i/. The /x/ 
dependence of the Binder cumulant used to determine the nature of the transition is very 
weak for the lattice sizes we use (8^ x 4, 12^ x 4 and 16'^ x 4). For this reason, our results can 
only be considered suggestive, and not definitive. Similar conclusions have been drawn by de 
Forcrand et al. from simulations at imaginary /x j^, 26 1. This disagrees with the early work 
of the Bielefeld-Swansea collaboration [21], who did claim to find such a critical endpoint. 
However, as indicated before, these simulations used the R algorithm which as de Forcrand 
and Philipsen and we discovered can lead to misleading results. 

De Forcrand et al. have recently introduced reweighting methods which enabled them 
to calculate the slope of the Binder cumulant directly, thus reducing the errors to a point 
where the sign is determined unambiguously This shows that the critical mass does 
indeed decrease with increasing /x, so that there is no critical endpoint associated with mc(0). 
However, their published results using this new method are all on 8^ x 4 lattices where finite 
size effects, such as the fact that the chiral condensate is not the true order parameter (in 
the renormalization group sense), are large. 

For this reason we have been investigating the use of such reweighting techniques for 
phase-quenched QCD on 12^ x 4 lattices. Larger lattices are less suited to such reweighting 
because the overlap between the ensembles of configurations at /x/ and /ij -|- for given 
/i/ and is smaller for larger lattices. The ratio of determinants is further from unity 
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for the larger lattices, and the fluctuations associated with the noisy estimator on a single 
configuration are also larger. Correlations in molecular-dynamics time are longer on the 
larger lattice. Our tests are promising and suggest that using finite rather than zero /i/ 
configurations for the reweighting are preferable. However, since reweighting is expensive, 
unless we can find a way to make better use of the fact that the ratio of fractional powers 
of Dirac operators for fij and fij + Sfij is better conditioned than either of the original 
operators, reweighting will be considerably more expensive than the cost of producing a 
single trajectory, so that it is unclear as yet whether reweighting will prove to be the most 
cost-effective method of getting definitive results on these larger lattices. 

One might ask whether our failure to find a critical endpoint disagrees with the work 
of Fodor and Katz [3]. They reported a critical endpoint at /i^ = 360(40) MeV and T = 
162(2) MeV, and hence fij = 240(27) MeV. Since our method breaks down for fij > m-,,: ^ 
140 MeV, their value is beyond the reach of our method. Hence our simulations do not show 
the absence of a critical endpoint, only the absence of a critical endpoint associated with 
the critical point at fi = fij = 0, for 3-fiavour QCD. 

We are now extending our simulations of 3-fiavour phase-quenched lattice QCD to enable 
a calculation of the equation-of-state of this theory outside the superfiuid region. This will 
enable comparison with full QCD. In addition, the phase-diagram of QCD at finite isospin 



chemical potential and its equation-of-state are o 



to new activity in the studies of these theories 321] 



interest in their own right. This has led 



All our simulations have been performed with the standard staggered action, with A^^^ = 4 
and are thus subject to large discretization errors. The Bielefeld-Swansea collaboration found 
that the critical mass at /i = decreased dramatically, when they changed their lattice action 
from the standard lattice action to a highly improved action, indicating that this mass is 



very sensitive to finite lattice-spacing errors 2l|, |22|. A less dramatic decrease in the critical 
mass has recently been reported by de Forcrand, Kim and Philipsen when they increased 
Nt from 4 to 6 with the standard staggered action 2^ ■ Hence we should consider repeating 
our simulations at larger Nt, improving the action we use, or both. 

In using a staggered action for 3 fiavours, we are ignoring the so-called 'rooting' contro- 
versy. People have questioned whether taking fractional powers of the fermion determinant 
to allow use of staggered fermions to simulate numbers of fermion fiavours which are not 
multiples of 4, defines a theory with a sensible continuum limit. We direct the reader to 
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33l . |3J] for recent arguments on both sides of this controversy. Even assuming that this 



controversy is resolved and indicates that 'rooted' staggered fermions are legitimate, this is 
only relevant to the case of zero fi. At ^ ^ 0, Golterman, Shamir and Svetitsky have pointed 
out that further ambiguities arise with regard to taking fractional powers of the phase of 
the determinant 35(]. We have avoided this difficulty by ignoring the phase, but eventually 



it will need to be faced. 
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